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We show that a reformulation of the ADM equations in general relativity, which has dramatically 
improved the stability properties of numerical implementations, has a direct analogue in classical 
electrodynamics. We numerically integrate both the original and the revised versions of Maxwell's 
equations, and show that their distinct numerical behavior reflects the properties found in linearized 
general relativity. Our results shed further light on the stability properties of general relativity, 
illustrate them in a very transparent context, and may provide a useful framework for further 
improvement of numerical schemes. 
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Motivated by the prospect of gravitational wave de- 
tections and the accompanying need for theoretical grav- 
itational wave templates, much effort has recently gone 
into the development of numerical relativity algorithms 
that are capable of modeling the most promising sources 
of gravitational radiation, in particular the inspiral and 
coalescence of binary black holes and neutron stars. In 
the past, progress has been hampered by numerical insta- 
bilities that arise in straight-forward implementations of 
the traditional Arnowitt-Deser-Misner (ADM, Jl)) 3 -f 1 
decomposition of Einstein's equations (e.g. ^). These 
instabilities have been associated with the mathematical 
structure of the ADM equations, and as a cure a number 
of hyperbolic formulations have been suggested (e.g. ^ 
as well as [D and references therein). Alternatively, Shi- 
bata and Nakamura Q and later Baumgarte and Shapiro 
suggested a reformulation of the ADM equations that 
has been demonstrated to dramatically improve the sta- 
bility properties of numerical implementations (e.g. [^). 
While the exact reason for this improvement is still some- 
what mysterious (see 1^, hereafter AABSS, and fl^), this 
new formulation, now often called the BSSN formulation, 
is quite widely used. 

In this Brief Report we show that the BSSN reformula- 
tion of the ADM equations has a direct analogue in clas- 
sical electrodynamics (E&M) . We numerically implement 
both the original and the revised versions of Maxwell's 
equations, and find that their distinct numerical prop- 
erties reflect those found in general relativity (GR). As 
suggested by AABSS, these properties can be identified 
with the propagation of constraint violating modes. We 
present our findings in the hope that, in addition to being 
a very transparent illustration of the stability properties 
of GR, they may prove useful for the future development 
and improvement of numerical algorithms. 

In a 3 -t- 1 decomposition of GR, Einstein's equations 
split into the two constraint equations 

R- K,jK'^ + = 2p (1) 

and 

D,K'^ - D,K = S, (2) 



and the two evolution equations 

dtlij = -2aK,j (3) 

and 

dtK.j = -D,D,a + a{R,j - 2KaK\ -f KK.^ - My ). 

(4) 

Here is the spatial metric, Kij the extrinsic curvature, 
K = ■y^^Kij its trace, a and /?' are the lapse function 
and the shift vector, and p, Si and Mij are matter source 
terms. The time derivative is defined as dt — dt ~ Cp, 
and Rij 

Rij = --jl'''' {lijM + IkLij - lki,il - lil,kj^ (5) 

Di and F*,. are the Ricci tensor, the covariant derivative, 
and the connection coefficients associated with 7^ . Fi- 
nally, R is the scalar curvature R = Y^Rij- Equations 
(|l|) through (^) are commonly refered to as the ADM 
equations 0]. 

The first term in the Ricci tensor (^ is an elliptic op- 
erator acting on the components of the spatial metric 
7ij. If the Ricci tensor contained only this term, the 
two evolution equations (||) and (Q) could be combined 
to form a wave equation for 7^- . This property is spoiled 
by the appearance of the three other second derivative 
terms in (^ , suggesting that these terms may be respon- 
sible for the appearance of instabilities in many straight- 
forward, three-dimensional implementations of the ADM 
equations. This problem can be avoided by either using a 
hyperbolic formulation of GR (e.g. ||]), or by eliminating 
the mixed second derivatives as in the BSSN formulation. 

In this formulation, the conformally related metric 7^ 
is defined as 7^ = e~'*'^7ij, where the conformal factor 
e"^ is chosen so that the determinant of •jij is unity, 7 = 
1. The conformal exponent (f) as well as the trace of 
the extrinsic curvature, K, are evolved as independent 
variables. Their evolution equations can be found from 
the traces of equations (||) and (||), while the trace-free 
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parts of those equations form evolution equations for 7^ 
and the trace-free part of the extrinsic curvature, Aij. 
The latter equation still contains the Ricci tensor Rij 
associated with 7^ , which contains all the mixed second 
derivatives of (||) . The crucial step is to realize that these 
second derivatives can be absorbed in a first derivative 
of the "conformal connection functions" 



(6) 



where the last equality holds because 7 = 1. Here and 
in the following all barred quantities are associated with 
7y . In terms of , the Ricci tensor can be written 

(compare Evidently, the only remaining second 

derivative term is the elliptic operator 7'™7ij,im, if the f * 
are considered as independent functions. For that pur- 
pose, an evolution equation is derived by permuting a 
time and space derivative in 



(8) 



The divergence of the extrinsic curvature can now be 
eliminated with the help of the momentum constraint 
(0), which yields the evolution equation 



(9) 



As in 0] we will refer to the ADM equations (|l|) through 
as System I, and to the new BSSN system of equations 
as System IF A complete listing and derivation of the 
latter can be found in 0. 

In an effort to better understand the improved numer- 
ical behavior of System II, AABSS linearized Systems I 
and II and identified their characteristic structure. They 
found that System I has constraint violating modes with 
a characteristic speed of zero. In System II, the charac- 
teristic speed of these modes changes to the speed of light. 
AABSS further demonstrated that in a non-linear model 
problem the existence of non-propagating, constraint vi- 
olating modes may lead to numerical instabilities, as en- 
countered in implementations of System I. Their analysis 
also demonstrated that the usage of the momentum con- 
straint in the derivation of is crucial for the charac- 
teristic speed for the constraint-violating mode to change 
to a non-zero value, and hence for the stability of the 
system. In the following we will show that very similar 
properties can be found in E&M. 



In terms of a vector potential Ai , Maxwell's equations 
can be written as the evolution equations 

dtA, = -E,-D,iP (10) 
dtE, = -D^D,A, + DW,A,~A7Tj,, (11) 



and the constraint equation 



(12) 



Here is the electrical field, pe the charge density, ji 
the flux, and ip the scalar gauge potential. Identifying 
Ai with 7y , Ei with Kij , and ^ with / ?', we see that 
the structure of the evolution equations ( p^ and ( [III ) is 
very similar to that of equations (||) and (|4|), and that 
the constraint equation (jlj) can be similarly identified 
with the momentum constraint (^. In analogy with GR, 
we refer to equations (|l^) through ( [T^ ) as System I. 

In further analogy, we will eliminate the mixed second 
derivative in (|ll|) by introducing a new variable F 



F = D,A' (13) 

(compare (||)), in terms of which (^ij) becomes 

dtE, = -DWjA, + D,r - 47rj,. (14) 

As in (^, the mixed second derivatives have been ab- 
sorbed in a first derivative of the new variable. An evo- 
lution equation for F can again be derived by permuting 
a time and space derivative in the definition of the new 
variable ( p^ ) 

dtT = DAA' = -D,E' - D'D,-<P = -Airp, ~ D'D,^,. 

(15) 

Here we have used the constraint (|l2| ), which, as in GR, 
will turn out to be crucial (see (|2l|) below). Equations 
(|lO|), (p^, (^5|) form the evolution equations of what we 
call System IF The definition (|l^) together with ( p^ ) 
form the constraint equations of System II. 

For vanishing scalar potential ■0 = 0, an analytical 
vaccuum solution to Maxwell's equations is given by the 
purely toriodal dipole field 



A-^ = ^sin6' 



2A 



ve 



(16) 



Here A is the amplitude, A parameterizes the size of the 
wavepackage, and u and v are the retarded and advanced 
time u = t + r and v = t — r. According to (|l0|), E^ can 
be found by taking a time derivative of (ITq) . Since p^ — 
and DiA^ — in this solution, ip ~ is consistent with 
both the Coulomb gauge and the Forentz gauge 



(17) 



(where the second equality applies for System II), if ap- 
propriate boundary and initial conditions are chosen. All 
results shown below were obtained with Forentz gauge. 
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FIG. 1: The integrated constraint violation ||C|| = (J C^dVY^^ for evolutions using Lorentz gauge with resolution 32^ and 64^ 
with outer boundaries at 4 (OB4) and 8 (OB8). 



As initial data for our dynamical simulations we adopt 
the analytical solution ([l^) at i = 



^"^ = 0, E'>' ^^ArsinOX^ 



(18) 



with ^ = A 1 and transformed into cartesian coordi- 
nates. 

We numerically implement System I and II following 
the algorithm of as closely as possible. In particular, 
we wrote the code in three spatial dimensions using carte- 
sian coordinates. We use an iterative Crank-Nicholson 



scheme |12 to update the evolution equations, and im- 
pose outgoing wave boundary conditions on E^, and 
r as in [0. We verified that the numerical solution con- 
verges to the analytical solution to second order as long 
as the solution is not affected by the outer boundaries 
(which are first order accurate). 

We compare the performance of System I and II by 
monitoring the constraint violation 



C = D,E' - Attp,. 



(19) 



In Fig. |l[ we show integrated values ||C|| = (/C^dV")^/^ 
for Systems I and II for two different gridsizes (32'^ and 
64'^) and two different locations of the outer boundary 
(at 4 (0B4) and 8 (0B8)). 

At early times, System I violates the constraint ( p^ ) 
to a lesser degree than System II. After about a light- 
crossing time, when the electro-magnetic wave has left 
the numerical grid, ||C|| settles down to a nearly constant 
value, which primarily depends on the grid-resolution. In 
System II, ||C|| is also largest after about a light-crossing 
time, but after that ||C|| decreases exponentially. As one 
might expect, the decay time of this exponential fall-off 
scales with the location of the outer boundary. 

To further compare these Systems, we compare snap- 
shots of the constraint violation C for the 32"^ 0B4 evolu- 
tion in Fig. ||. With identical initial data, both Systems 
have identical values for C at i = 0. At an intermediate 



time t — 9.375, C is larger in System II than in System 
I, as one expects from Fig. ^ At a much later time 
{t = 200), however. System I has settled down to a con- 
stant shape, while C in System II has almost completely 
dissipated. 

These numerical results demonstrate that, as in GR, 
the constraint violation C behaves very differently in the 
two systems. In E&M, this different behavior can be 
understood very easily from an analytic argument. 

For System I, it is easy to show that a time derivative 
of the constraint violation (|l^) vanishes 

dtC = DAE' - AiidtPe - -^AD^f + dtPe) = 0, (20) 

where we have used the continuity equation Dif" +dtPe = 
0. This explains why at late times the profile of C remains 
unchanged in System I. This property is the analogue of 
AABSS's finding that the linearized ADM equations have 
non-propagating, constraint violating modes. 

For System II, on the other hand, it can be shown that 
C satisfies a wave equation 



= dtD'dtE,~iTTdfpe 

= dtD\-DjD'JA, + AT - 4^ji) - And^pe 

= -D\D,D^dtA, - D.dtT) - A7rdt{D\n + dtPe) 

= D' (DjD^{E, + D,4,) - A(4^Pe + DWj^)^ 

= DjD'{D'E, - Aixpe) = DjD'C, (21) 



which explains why the constraint violations propagate 
away in System II. This is the analogue of AABSS's find- 
ing that in the relativistic System II (the BSSN equa- 
tions) the constraint violating modes propagate with the 
speed of light. Moreover, wc now realize why the usage 
of the constraint ( p^ ) in (^ was crucial - without hav- 
ing made this substitution the terms on the last line of 
( pl| ) would have cancelled, dfC = 0, leading to a non- 
propagating constraint violation as in System I. The ad- 
dition of constraint equations to the evolution equations 



FIG. 2: C in the z = .0625 plane (the first of our grid points) for System I (left column) and System II (right column) for a 
32^ OB4 evolution in Lorentz gauge. We show results for t — Q (top row), t — 9.375 (middle row) and t = 100 (bottom row). 
C settles down to a constant profile in System I, but dissipates away in System II. 



has been found to be crucial in several other formulations 
of Einstein's equations as well (e.g. p^). 

To summarize, we see that the numerical stability 
properties of two formulations of GR are beautifully re- 
flected in similar formulations of E&M. Maxwell's equa- 
tions therefore provide a very transparent framework for 
analyzing these properties, which may be useful for fu- 
ture algorithm development. For example, Fig. |^ shows 
that the outer boundaries produce much more noise in 
System I than in System II, which points to an incon- 
sistency between the treatment of the interior equations 



and the boundary. Similar problems are likely to occur 
in GR, but might be easier to analyze in the simpler 
framework of E&M. 
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